#!/usr/bin/perl
# 
# Exercise 9.5
# 
# Include a graphic display of the cut sites in the restriction map by printing the 
# sequence and labeling the recognition sites with the enzyme name. Can you make a 
# map that handles multiple restriction enzymes? How can you handle overlapping 
# restriction sites?
#
# Answer to Exercise 9.5 - Part 1

use strict;
use warnings;


# This is a tricky exercise.  The beginning of it is fairly simple.  You just need
# to define a "commentary" array of lines, and place the restriction enzyme names at the
# corresponding locations: then print the lines of the sequence and the lines of the
# array interleaved.  So let's start with that.  We'll use Exercise 9.3, but add a
# subroutine "graphicrestrictionmap1" to handle the output.  This is a bit simple-minded
# - for instance, what if the enzyme name gets split between two lines?  But it
# does the basic job.
# 

# Declare and initialize variables
my %rebase_hash = ();

unless(dbmopen(%rebase_hash, 'DBNAME', 0644)) {

	print "Cannot open DBM file DBNAME with mode 0644\n";
}

my @file_data = ();
my $query = '';
my $dna = '';
my $recognition_site = '';
my $regexp = '';
my @locations = ();

# If there is a command-line argument, assume it's DNA
if(@ARGV) {
    $dna = $ARGV[0];

# Otherwise, prompt for a FASTA file to open
}else{

    print "Input a FASTA filename: ";
    my $filename = <STDIN>;
    chomp $filename;

    unless(-e $filename) {
    	print "$filename does not exist!\n";
	exit;
    }
    # Read in the file 
    @file_data = get_file_data($filename);

    # Extract the DNA sequence data from the contents of the FASTA file
    $dna = extract_sequence_from_fasta_data(@file_data);
}

# Get the REBASE data into a hash, from file "bionet"
%rebase_hash = parseREBASE('bionet.110');

# Prompt user for restriction enzyme names, create restriction map
do {
    print "Search for what restriction site (or quit)?: ";
    
    $query = <STDIN>;

    chomp $query;

    # Exit if empty query
    if ($query =~ /^\s*$/ ) {

        exit;
    }

    # Perform the search in the DNA sequence
    if ( exists $rebase_hash{$query} ) {

        ($recognition_site, $regexp) = split ( " ", $rebase_hash{$query});

        # Create the restriction map
        @locations = match_positions($regexp, $dna);

        # Report the restriction map to the user
        if (@locations) {
            print "Searching for $query $recognition_site $regexp\n";
            print "A restriction site for $query at locations:\n";
            print join(" ", @locations), "\n";

	    ##################################################
	    #
	    # Here is the modification to Exercise 9.3
	    #
	    ##################################################
 
	    my @output = formatrestrictionmap1($dna, $query, @locations);

	    print @output;

	    ##################################################
	    #
	    # Here ends the modification
	    #
	    ##################################################

        } else {
            print "A restriction enzyme $query is not in the DNA:\n";
        }
    }
    print "\n";

} until ( $query =~ /quit/ );

exit;

################################################################################
# Subroutines
################################################################################

################################################################################
#
# Subroutine match_positions
#
# Find locations of a match of a regular expression in a string
#
# 
# return an array of positions where the regular expression
#  appears in the string
#

sub match_positions {

    my($regexp, $sequence) = @_;

    use strict;

    #
    # Declare variables
    #

    my @positions = (  );

    #
    # Determine positions of regular expression matches
    #
    
    while ( $sequence =~ /$regexp/ig ) {

        push ( @positions, pos($sequence) - length($&) + 1);
    }

    return @positions;
}

################################################################################
#
# Subroutine formatrestrictionmap1
#
# Interleave a sequence with annotation lines showing the locations of an enzyme
#

sub formatrestrictionmap1 {

    my($seq, $enz, @pos) = @_;

    # In Chapter 8 we had a subroutine print_sequence that formatted and printed
    # sequence data.  Now maybe we wish we'd made it more modular, and that we
    # had a subroutine that only formatted a sequence, returning an array to be
    # printed.  Here we use our new subroutine "format_sequence" - the definition
    # appears later in the code.

    my $line_length = 50;
    my @formatted_sequence = format_sequence($seq, $line_length);

    # Now we have to take the @pos of enzyme positions, and create a new array
    # that contains an array of spaces exactly as long as the array
    # @formatted_sequence, with the enzyme name inserted at the right locations.
    # Here's a simple way to do it: let's create a scalar of spaces the same length
    # as the sequence; insert the enzyme name at the appropriate places; then split
    # the scalar into an array using our format_sequence subroutine.  Here goes:

    # Make an "blank" annotation line of the correct length
    my $annotation = ' ' x length($seq);

    # Put the labels for the enzyme name at the correct positions in the annotation
    foreach my $location (@pos) {
        substr($annotation, $location-1, length($enz)) = $enz;
    }

    # Make a formatted array of the annotation
    my @formatted_annotation = format_sequence($annotation, $line_length);

    # Merge the formatted annotation with the formatted sequence

    my @output = ();

    while(@formatted_annotation) {
        push(@output, shift @formatted_annotation);
        push(@output, shift @formatted_sequence);
	# separate the lines
	push(@output,"\n");
    }

    # Return the merged annotation and sequence
    return @output;
}
    

sub format_sequence {

    my($seq, $line_length) = @_;

    my(@output) = ();

    # Split sequence into lines of $line_length
    for ( my $pos = 0 ; $pos < length($seq) ; $pos += $line_length ) {
        push(@output, substr($seq, $pos, $line_length) . "\n");
    }

    return @output;
}

